clear; close all; clc;

%% 参数设置
alpha = 5;
sigma = 0.8;
mu = 0.1;
epsilon = 0.1;   % φ 更新系数
k = 0.2;         % 固定耦合强度k

N_total = 2000;  % 总迭代次数
N_cut = 1000;     % 舍去的暂态迭代次数

% φ0取值范围
phi0_min = -2;    
phi0_max = 6;     
N_phi0 = 1000;    
phi0_values = linspace(phi0_min, phi0_max, N_phi0);

% 创建空数组用于存储极值部分的数据
Phi0_1 = [];
X1_record = [];
Phi0_2 = [];
X2_record = [];

for idx = 1:length(phi0_values)
    phi = phi0_values(idx);
    [X1_save, X2_save, ~] = Rulkov(alpha, sigma, mu, k, phi, N_total, epsilon, N_cut);
    
    % 取稳态数据
    X1_steady = X1_save(N_cut+1:end);
    X2_steady = X2_save(N_cut+1:end);
    
    % 提取X1极值（峰谷值）
    [peaks1, ~] = findpeaks(X1_steady, 'MinPeakProminence', 0.05);
    [troughs1, ~] = findpeaks(-X1_steady, 'MinPeakProminence', 0.05);
    troughs1 = -troughs1;  % 还原谷值
    
    % 确保 peaks1 和 troughs1 为列向量，并合并
    extrema1 = [peaks1(:); troughs1(:)];
    
    % 提取X2极值
    [peaks2, ~] = findpeaks(X2_steady, 'MinPeakProminence', 0.05);
    [troughs2, ~] = findpeaks(-X2_steady, 'MinPeakProminence', 0.05);
    troughs2 = -troughs2;  % 还原谷值
    
    % 确保 peaks2 和 troughs2 为列向量，并合并
    extrema2 = [peaks2(:); troughs2(:)];
    
    % 保存数据
    Phi0_1 = [Phi0_1; phi*ones(length(extrema1), 1)];
    X1_record = [X1_record; extrema1];
    Phi0_2 = [Phi0_2; phi*ones(length(extrema2), 1)];
    X2_record = [X2_record; extrema2];
end

%% 绘制分岔图（只绘制极值部分）
figure;
hold on;
scatter(Phi0_1, X1_record, 1, 'b', 'filled');  % X1 用蓝色标记
scatter(Phi0_2, X2_record, 1, 'r', 'filled');  % X2 用红色标记
xlabel('phi');
ylabel('X1, X2');
title('分岔图');
legend('X1', 'X2', 'Location', 'best');
hold off;